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Simulating Bosonic Baths with Error Bars 
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We derive rigorous truncation-error bounds for the spin-boson model and its generalizations to arbitrary 
quantum systems interacting with bosonic baths. For the numerical simulation of such haths the truncation of 
both, the number of modes and the local Flilbert-space dimensions is necessary. We derive super-exponential 
Lieb-Robinson-type bounds on the error when restricting the bath to hnitely-many modes and show how the 
error introduced hy truncating the local Hilbert spaces may be efficiently monitored numerically. In this way 
we give error bounds for approximating the infinite system by a finite-dimensional one. As a consequence, 
numerical simulations such as the time-evolving density with orthogonal polynomials algorithm (TEDOPA) 
now allow for the fully certified treatment of the system-environment interaction. 


Introduction - Ideal quantum systems may be considered 
closed, undergoing textbook unitary evolution. In any realis¬ 
tic experimental setup however a quantum system is open, that 
is, it interacts with an environment composed of those degrees 
of freedom that are not under the control of the experimenter. 
Hence the numerical and analytical description of the dynam¬ 
ics of a quantum system in interaction with its environment 
is of fundamental importance in quantum physics. The pre¬ 
cise nature and composition of the system-environment inter¬ 
action is generally not known, but for a wide range of systems 
encountered in physics, chemistry, and biology, it is common 
to model the environment as a continuum of harmonic oscil¬ 
lators, which interact linearly with the system. This results 
in the paradigmatic spin-boson model that captures many as¬ 
pects of the system-environment interaction [1]. The spin- 
boson model is exactly solvable only in the rarest of special 
cases and one is therefore compelled to employ a variety of 
approximations and numerical descriptions in order to ob¬ 
tain the reduced dynamics of the quantum system in question. 
Notable examples include those cases in which the environ¬ 
ment possesses a correlation time that is much shorter than 
the system dynamics and the system-environment interaction 
is weak. Under these assumptions it is then well-justified and 
customary to resort to the so-called Markov approximation 
which permits the derivation of completely positive and linear 
differential equations, the Lindblad equation, for the quantum 
system alone [2]. 

However, settings of considerable practical importance may 
violate either or both of these assumptions and require a 
more sophisticated treatment. The recently emerging inter¬ 
est in quantum effects in biological systems provides a case 
in point [3]. For instance, in typical pigment-protein com¬ 
plexes the dynamical time-scales of the vibrational environ¬ 
ment can be comparable or even slower than the quantum me¬ 
chanical excitation energy transfer dynamics. Moreover, in 
the limit of slow bath dynamics, perturbative treatments of 
the coupling between system and environment cannot be used 
even if the system-bath coupling is intrinsically weak. Con¬ 
sequently, steps have been taken towards the development of 
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FIG. I. A system coupled to a bosonic bath. Red lines indicate the 
truncations: The spatial truncation to a chain of finite length L and 
the truncation of the local Hilhert space dimensions to rrii. 


non-perturbative and non-Markovian approaches for the de¬ 
scription of the quantum system-environment interaction (see 
[3, 4] for overviews of recent developments). However, the 
majority of these approaches have in common that they exploit 
approximations that are not well controlled in the sense that no 
rigorous error bounds on the simulation results are available. 
Hence these methods are not certified. 

The time evolving density with orthogonal polynomials al¬ 
gorithm (TEDOPA) for the spin-boson model presents a no¬ 
table exception, as will be demonstrated in the present work. 
It makes use of an exact transformation of the standard rep¬ 
resentation of the spin-boson model onto a spin interacting 
with the first site of a semi-infinite nearest-neighbor coupled 
chain [5-9] which renders the system particularly amenable 
to time-adaptive density matrix renormalisation group (t- 
DMRG) simulations. The structure of the resulting system is 
such that excitations tend to propagate along the chain away 
from the system towards infinity leading to irreversible sys¬ 
tem dynamics for long times. This approach has been used 
with success in the simulation of a number of highly non- 
Markovian system-environment interactions [6, 10, 11]. 

The errors that accumulate in the t-DMRG simulation can 
be bounded rigorously. Nevertheless, the numerical TEDOPA 
simulation employs two as yet uncertified assumptions: (i) 
the semi-infinite chain needs to be truncated to a finite length 
and (ii) the local dimension associated with each harmonic 
oscillator of the chain the needs to be truncated to a finite di¬ 
mensional Hilbert space, see Pig. 1. The errors that are intro¬ 
duced in this manner are usually estimated by increasing both 
the chain length and Hilbert space cut-off until the change in 
the result drops below a predefined threshold. However, in 
practice this somewhat inelegant approach can become highly 
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challenging numerically, and can lead to erroneous numeri¬ 
cal predictions [12]. A more rigorous approach is therefore 
desirable. 

Here we employ techniques that lead to Lieb-Robinson 
type bounds to achieve this goal by deriving bounds for the 
errors arising from approximations (i) and (ii). As the errors 
arising in each step of the t-DMRG integration can also be 
bounded we arrive at a method that possesses rigorous error 
bounds on the results that it delivers. This extends signifi¬ 
cantly existing recent results in the literature that apply to the 
finite dimensional setting of spin systems [13] and therefore 
allows the fully certified treatment of the system-environment 
interaction for both, harmonic oscillator as well as spin envi¬ 
ronments. 

The system under consideration - We will consider the 
Hamiltonian of an arbitrary system Hs coupled via y to a 
bosonic bath described by Hb so that the total Hamiltonian 
reads 


H = Hs + V + Hb. ( 1 ) 

For simplicity and to directly connect to the TEDOPA ap¬ 
proach [6, 7, 10, 11], we assume that Hb describes a one¬ 
dimensional nearest-neighbour Hamiltonian (the higher di¬ 
mensional case with more general couplings will be published 
elsewhere [14]) and takes the form 

1 °° 

Hb — ^ ^ ^ piPijPj '), ( 2 ) 

ij=0 


conditions is satisfied). Let c be such that < c. 

Then 


4||OP||/,||/c 




-f 1) 
(L + 1)! 


, ( 5 ) 


where C = \\Pl\\\Xl-i,l \/+ \Pl-i,l\/ c and 


7o 


Ixx ^xp 
Ipx Ipp 


[lab]i,j = tr[ai6jPo], 


( 6 ) 


collects the two-point bath correlations in the initial state. If 
P (X % we may replace L by 2L in Eq. (5). 

If the initial 2-point correlation functions (the matrix elements 
of 7 o) are unbounded, then one can still achieve bounds, see 
the appendix for details. The r.h.s. of Eq. (5) describes the 
Lieb-Robinson-type light cone [15]. Outside the light cone, 
so for r := ect < L, one finds super-exponential decay in L: 
{ct)^e^*/L\ < This makes rigorous the physi¬ 

cal intuition that for all finite times only a chain of finite length 
is required to simulate the dynamics of local observables to 
within a prescribed precision. Our bound applies to any sys¬ 
tem Hamiltonian, unbounded or otherwise, and depends only 
linearly on the operator norm of the system coupling || /i||. The 
proof relies on Lieb-Robinson bounds for harmonic systems 
[16-18] (see also Ref. [19]) and may be found in the appendix. 
Before stating our second main result, we discuss the above 
bound in the light of the generalized spin-boson model. 

Generalised spin-boson model - In this section we will in¬ 
vestigate Hamiltonians of the form 


where we assume that only nearest-neighbours are coupled, 
Xij = Pij = 0 for \i — j\ > 1, and we let w.l.o.g. 
Xij = Xji G R, Pij = Pj^i G R. We consider system- 
bath couplings of the form V = h® xq (see the appendix for 
systems coupled to several baths), where h acts on the system 
and we assume that it is bounded in operator norm, || /i|| < oo. 
The system with Hamiltonian Hs has no restrictions, it can 
correspond to any system—^bosons, fermions, and/or spins, all 
in arbitrary dimensions. 

Spatial truncation of the bath - Eor bounded system ob¬ 
servables O, ||0|| < oo. We are interested in the quantity 

A(/,L) = |tr[6e-'^‘poe‘^‘] - tr[6e-‘*^‘poe*"^*] |, (3) 

i.e., the error introduced when, instead of simulating the full 
Hamiltonian H , we simulate the time evolution of system ob¬ 
servables O with the truncated bath Hamiltonian 

1^-1^ 

Hb = ^ ^ ' i^iXijXj -f piPijPj) (4) 

i,j=0 

and corresponding total Hamiltonian Hl = Hs V -\- Hg. 
Our first main result is the following. 

Theorem 1 Let H and Hb be as above. Let X,P > 0 or 
X = P (see the appendix for a bound when neither of these 


H = Hs + 


dkg{k)alak-\-As /d. h{k){al-\-ak). (7) 


This describes a quantum system with Hamiltonian Hs inter¬ 
action with a bath of bosons; it is described in more detail in 
terms of second quantised operators in [20]. This model has 
received renewed interest in recent years due to its importance 
in the theoretical study of quantum effects in biology (see [3] 
for a review). An important quantity that describes the bath 
and its coupling to the system is the spectral density, which, 
for invertible g, is defined as 


J{u}) = irh'^ (g ^(w)) 


dg 

du! 


( 8 ) 


with g~^ the inverse of g. The smallest closed interval con¬ 
taining the support of g~^ is denoted The case 

bOmin = 0 is Called massless where as Wmm > 0 is known as 
massive. 

Building on the work of [5, 7, 8], it was shown using the 
theory of orthogonal polynomials in [9] that Eq. (7) can be 
written in the form of Eqs. (2,1) and that there are two ways 
to do this. Both choices 


h = poAs, X = P, (9) 

and 

h = P = OJrnax'^^ (10) 
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FIG. 2. Fock space truncation error (Eq. (17)) for the particle mapping and power-law spectral densities as in Eq. (19) with A/uJc = 1, 
a. — 0.8, s = 3 for initial state Po = <8 Q%, Q% = |t)(tl and q% the vacuum. We truncate each local Hilhert space at the same value 

rui = m and L has the values 3 to 6, but are indistinguishable (e.g. the difference between the L — 6 and L = 3 curve at the point denoted 
by a red square is 4.95 x 10“®). Lines are guides to the eye. The log-log plot on the left suggest algebraic increase in time and the plot on the 
right suggests better than exponential decrease with m. 


with appropriate X (given in terms of the spectral density in 
the appendix) are equivalent to Eq. (7). Here, 

ftp = - / dw J(w), ^— [ dw J(v/w) (11) 

and one finds ||X|| = ||P|| = ojmax for both cases and X > 0 
iff ijJmin > 0. Due to the form of their elementary excitations, 
the mappings leading to couplings as in Eqs. (9) and (10) were 
named particle mapping and phonon mapping, respectively, 
and we will adopt this denomination here. Crucially, in both 
cases, X couples nearest-neighbours only such that the bound 
in Eq. (5) is readily applicable to the particle and the massive 
phonon case, setting c = ujmax for both (similar results hold 
for the massless case, see appendix for full details). Eor the 
particle mapping, we hnd C < 2 and for the phonon mapping 
C < 1 such that, up to the constants /Tq/j, we obtain the same 
behaviour of the bound in both cases but replacing L by 2L 
in the massive phonon case. Hence, for the phonon mapping 
with a chain of only half the length, one has approximately 
the same chain truncation error as for the particle mapping. 

If the maximum frequency of the bath ujmax = oo, the 
chain coefficients are unbounded [9] and our bounds diverge. 
This divergence is not surprising in light of the observation 
that certain one-dimensional inhnite harmonic lattice mod¬ 
els with nearest neighbour interactions and unbounded coeffi¬ 
cients have been proven not to have a light cone bound [22]. It 
is noteworthy, that similar results can be derived for the case 
of a fermionic bath, since the chain mapping is still valid and 
Lieb-Robinson bounds for fermions are well-known [23]. 

Truncating local Hilbert spaces - We now consider the er¬ 
ror introduced when the local Hilbert space dimensions of the 
harmonic oscillators making up the bath are truncated. To this 


end, we define the projector 

m 

Hm. — ^rriQ '^m — ^ ^ 1^) (^17 (1^) 

n—0 

where acts on the i’th site of the bath and truncates the 
local Hilbert space according to 1^- Eor bounded observables 
acting on the system O, ||0|| < oo, we consider 

Xm{t)= |tr[(9e“**^poe**^]-tr[6e““^™poe“*”*]|, (13) 
i.e., the error introduced by evolving the system according to 

Hm = tmHlm (14) 

instead of H. Here, H is as in Eq. (4) and we omit the index L 
for notational clarity. The truncated Hamiltonian reads Hm = 
Hs + H]^ -\-h® tmXo'^m, where 

1^-1 

~ c\ ^ ^ Hm “F Pi,j'^mPi'Pj^m\ • (15) 

i,j=0 

In the appendix we show that 

< tr[(l - l^)po] + 2 / dx ^em{x), (16) 

4||0||^ Jo 

where 

em{x) = tr[/i^e““'^B (17) 

with 

X{x) = 

( 18 ) 

0m(a;) = 
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Crucially, under the assumption that the system Hilbert space 
is finite dimensional, this error may be computed numeri¬ 
cally as it involves only observables acting on the truncated 
Hilbert space is a linear combination of the 

Xi and Pi) and which are of a form amenable to t-DMRG sim¬ 
ulations (see the appendix for details). For all finite times, 
= 0 and we study its behaviour in m at the 
hand of numerical examples below. If the bath initially con¬ 
tains only a finite number of particles, tr[(l — lm)po] van¬ 
ishes for appropriate m. Such states include the vacuum state 
which is also the zero temperature thermal state for the par¬ 
ticle mapping. For higher temperature thermal states of the 
bath, tr[(l — l^)po] vanishes exponentially for large {rrii}. 
The total error induced on the expectation value of O due to (i) 
truncating the chain to finite length and (ii) the truncation of 
the local dimensions is bounded by the sum of the two individ¬ 
ual error bounds: A(f, L) + A^(f). This rigorously bounds 
the error of approximating an infinite-dimensional bath of 
bosons by a chain of length L made up of finite-dimensional 
subsystems with nearest neighbour interactions. If in addition 
we assume the system with Hamiltonian Hs to be a spin sys¬ 
tem, then the Hamiltonian is in the class which, as [13] shows, 
can be simulated with resources polynomial in L and error e, 
and exponential in |f|. 

Numerical example - As an example, we consider the spin- 
boson model with power-law spectral density, 

J(a;) = 7ra 0(1 — w/wc), (19) 

where 0 is the Heaviside step function. This model has been 
extensively probed numerically, and there has been contro¬ 
versy over the accuracy of numerically derived critical expo¬ 
nents. One of the issues with the results was the inability to 
verify the local Fock space truncation errors [12, 24]. The sys¬ 
tem Hamiltonian and interaction part are Hs = — A(Ta;/2 and 
= cr2/2. The dissipation is known as Ohmic for s = 1 


and super ohmic for s > 1. This can be written in the chain 
representation using Eq. (9) (see the appendix for details). In 
Fig. 2, the bound for the particle mapping is plotted for the 
super-ohmic case and various L and m. Constants used for 
the simulation (see figure caption) are taken from the literature 
[21]. The initial state of the bath corresponds to the zero tem¬ 
perature thermal state. We probe the same initial state for the 
case of ohmic dissipation and achieve qualitatively the same 
results (see appendix). Furthermore, we test the bound for a 
squeezed vacuum state of the bath, which is a highly popu¬ 
lated state (see appendix). 

Conclusion - The detailed simulation of the interaction of a 
quantum system with a structured environments composed of 
harmonic oscillators has applications in a wide variety of sci¬ 
entific fields. The multitude of proposed algorithms to tackle 
this problem numerically lacked a method that delivers a sim¬ 
ulation result with a rigorous error bound associated with it. 
In this work we derived error bounds that demonstrate that 
the recently developed TEDOPA can provide such a method. 
More specifically, obtaining Lieb-Robinson type expressions 
we provide complete error bounds on the simulation of ob¬ 
servables of quantum systems coupled to a bosonic baths with 
infinitely many degrees of freedom such as the spin-boson 
model. This includes the errors incurred due to the truncation 
of the local Hilbert-spaces of the harmonic oscillators and due 
to the truncation of the length of the harmonic chain repre¬ 
senting the environment. In this manner we provide a fully 
rigorous upper bound on the error for the numerical simula¬ 
tion of a spin-boson model and its generalisation to multiple 
baths and more general systems. 
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Appendix A: Spatial truncation of the bath 

We consider an arbitrary (not necessarily finite-dimensional) system described by a Hamiltonian Hs and a bosonic bath 
described by 


Hb 


CXJ 

ij=0 


(Al) 


where Xi and pi are canonical position and momentum operators with the usual commutation relations Xj] = [pi^Pj] = 0 
and [xi,pj\ = i()ij- As we are allowing the bath to consist of infinitely-many modes (infinitely-many lattice sites), we assume 
throughout that the domain of the Hamiltonian is well-defined. W.l.o.g., we let Xij = Xj^i S IR and Pij = Pj^i G R. We 
assume that they couple only nearest neighbours: Xij = Pij = 0 for \i-j\ > 1 - 
We suppose that system and bath are coupled according to 


V = hxo, 


(A2) 


where h acts on the system and we assume ||/i|| < oo. Thus, our total Hamiltonian reads 

H = Hs + V + Hb, (A3) 

where for compactness of notation, we have neglected tensor products with the identity. We are interested in the error introduced 
for the time-evolution of bounded observables O (assuming ||0|| < oo) acting on the system when, instead of simulating the 
full Hamiltonian H, we take only finitely many lattice sites of the bath into account. Namely those that are closest to the site 0, 
truncating the bath Hamiltonian according to 

1 1 » 

^B 2 ^ ^ “b Pi^i,jPj^ 2 ^ “b Pi{^L')i,jPj'\ i (A4) 

ij=0 i,j 

where Xl, Pl are the principle submatrices of X, P corresponding to the non-truncated modes. The truncated chain hence 
consists of L modes. Denoting the initial state of the whole system by qq and the total truncated Hamiltonian hy Hl = 
Hs + V + Hg, we set out to bound the difference 

A(f,L) = |tr[6e-‘^*poe‘^*] - tr[6e-'^^*poe‘^^‘] |. (A5) 

We will prove the following theorem. 


Theorem 2 (Spatial truncation of the hath) Let H, Hg as above, c, d such that < c and max{|| A||, ||P||} < d. 

Then 


A2(f,L)<4||Of^| 


|Pl_i,l|\ (cf) 


L+l 


(L + 1)! 


(e^‘ + l)(||7o||'/^ + 


- 1 




(A 6 ) 


If X, P > 0 or X = P, we may take c —>■ 0 such that we recover the theorem in the main text. If P oc 1, we may replace 

/ ^\2L + 1 

(2L+1)! - 


\ ixp IPP J 

collects the two-point bath correlations in the initial state of the whole system. Note that || Ai|| < ||A|| and ||Pl|| < ||P||- 

One can allow for the two-point correlations collected in 70 to diverge and still get a bound on A(f, L), see Section Ale. 
Often, one encounters systems interacting with multiple baths. We generalize to this setting in Section Aid. 
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1. Proof 


Denote 

U(t) = UL{t) = (A8) 

Then for system operators O 

tr[6e-'*^Poe“^] = tT[de-'*^^-^^U{t)goU^ = tr[e'*^^de-'*^^U{t)goU^t)] (A9) 

and similarly 

tr[C>e-'*^^Poe“^^] = de-'^^^UL{t)goUl{t)]. (AlO) 

Hence, 

tr[6e-'*^Poe'‘^] - tr[6e-*‘^"= tr[e‘‘^^6e-**^« (17(f)- 171(f)] + Wit) - 17i(f)]pol7l(f))] ■ (AH) 

Using the Cauchy-Schwarz inequality |tr[ABpo]P < ti'[AA^po]tr[^^^Po] ^ II Aptr[i?^i3po] the triangle inequality, we 
find 


A(f,L) = |tr[6e - tr[6e < 2||6||A7tr[[17t(f) - 17l(f)][(7(f) - 17L(f)]eo], (^^2) 


where 


and 


where 


tv[[U\t)-Ulit)][U{t)-ULit)]go] =-2^J dx^tr[U^ix)ULix)go] 


-i—UUx)UL{x) = 17Va;)e“^®e“'^« (V - ^ 

dx 


Y'^xHbYxH^Y^-IxH^YxHb _\r — I 


— V = —ih / dye 

Jo 




JvHb 


(A13) 


(A14) 


(A15) 


Let us summarize the bound so far; 


< fdx [dy |tr[i7l'(a;)e“^®e“'f^«e-'^'^^li[(i7B - |. (A16) 

8||0|p Jo Jo 


We now proceed to bound the commutator and come back to Eq. (A16) after Eq. (A18). We have 


L~1 


oo oo 


Hb — Hg = - [xiXijXj + piPi jPjj + - [xiXijXj + PiPijPj^ 

i—Lj—0 i—0j=L 

such that, as only nearest neighbours are coupled and X and P are symmetric. 


(A17) 






PB-i,e^^^-xoe-'^y^^ 


Pl-i,lPl 


=: C^^l_MXL-i,LiL + CI%.MPl-i,lPl- 

Let us now come back to Eq. (A16). Inserting the above expression, we see that we need to bound terms of the form 

Fr{x,y) = |_ 


(A18) 


(A19) 


with r = x,p. Writing ri(f) = g(x) = goe'^^’^e inserting the definition of 

U(x), and using [775,77 b] = [775,77^] = [775,7^(7)] = 0, this reads 


Frix,y) = |tr[e“(^^-^(e-“^^e“^iie-“^®e-“^^ri(a;-y)p(;r)]| < ||/i|| Jtr[f2 (x - y)p(x)] 


(A20) 













8 


where we used |tr[A_Bg(a;)]p < j] A|ptr[i?^Bg(a;)] to obtain the second line. Inserting Eqs. (A18,A20) into Eq. (A16), we hence 
have 


< \xl-i l 

8\\o\mh\\ - 


dx 




yhx{x,y) + \Pl-i,l 



yhp{x,y), 


(A21) 


where we denoted 


lr{x,y) = tr[fl{y)Q{x)\. (A22) 

To keep track of the case P cx 1, we let Pi j = 0 for \i — j\ > R with i? = 0,1. By Eq. (56) in Ref [16] and as ||AiiP£,|| = 
\\PlXl\\ [25], we have 


^ I |2n+l 

IQ,l-i(y)l < E 

n=0 ^ ' 

L<(n+l)(l+i?) 

(A23) 

n%-iiy)\< E 

L<ra(l + fl) + l 

Bounding the second moments 7r-(a;, y) in the following section, we return to Eq. (A21) in Section A 1 b to complete the proof. 


a. Second moments 


Recalling that q{x) = v)^ 

-i^tr[fi(y)p(a:)] = tr[fi(j/)[p(a;), 

= tr [p(a;)e“^^ , f i (y)]] 

= 2 tr [P"*- H (y)]fi (y)]. 


(A24) 


Now, 


fk{y) = = E cl%y)xi + E Ck,i(.y)pu 


(A25) 


where r = x^y and 


Cxxiy) Cxp{y) 
Cpx{y) Cpp{y) 


xpyyj \ ifs = AT 0 P, a = 


0 -1 
1 0 


(A26) 


Hence, 


L-l 


-i£tr[fi(y)p(x)] =2tr[p(x)e-^®/re--"-fi(y)] E 0 c2^,(y)w] 


z=o 

L-l 


2itr[p(x)e-^4e--*-fi(y)] E {cZ{y)d^o:ii^) - cZ{y)dZ{x)) , 


1=0 


(A27) 


where 


dxrix) drn(x) \ ttT ^ / 0 —'^L 

xx\ } xp\ I \ aLtisX^ Hb=Xl®Pl, CTl = ' ^ 


dpx(^x^ dpp{x^ 


11 0 


(A28) 


We find 


q-^HbV^ 


U 0 

0 0 I -aBH-x ( U 0 0 0 

0 1l r V 0 0 li 0 
0 0 


^ Cxpiy) i ) - Cxx{y) 


dip{x) 0 \ \ 

0 0 / 0 0 

/ X ( dLJx) 0 \ , , / dl.Jx) 0 , 

c™(y) ( 0 0 J ~V 0 0 j */ 


(A29) 
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the operator norm (and therefore the absolute value of all entries) of which is upper bounded by ||e^®°’^|| ||. Therefore, 

employing \a[ABQ{x)\\^ < |jAptr[B'l'Bp(x)], 

^ tr[fi(y)p(a:)]| < 2||e^«‘"^|| tr[p(a;)e“^®/ie-“^sfL(t/)] < 2||e^^"^|| 


dx 

which implies 


lr{x,y) < ^JtY[fl{y)go\ + ||/i ||f dz ||e 




tT[g{x)rl{y)], 

(A30) 

(A36) 


From Eq. (A25), we have, denoting ['yxxiy)]i,j = ir[^iiy)xjiy)go], ['yxpiy)]i,j = tAxiiy)Pjiy)Qo], hpx{y)]t,j = 
tr[My)xjiy)go], ['lppiy)]i,j = tr[pi{y)pj{y)go], 




(A37) 


i.e., all entries are bounded from above by || 7 o|| in particular [jrr{y)]L,L such that 


7.(^,y)<ll7oir/^||e^-‘^^|| + ||^||||e^-‘^^|| / dz (A38) 

Jo 

If Hb > 0, we may use the Williamson normal form to write Hbct = S^{D 0 D)Say = S^{D 0 D)a{S*)~^y, where 
{D 0 D)a is real skew-symmetric, i.e., its eigenvalues are purely imaginary. Hence, ||e^^'^^|| = = 1. We also 

have Ije^^'^^ll = 1 if X = P as then aHs is real skew-symmetric. Hence, c' an upper bound to max{||X||, ||P||}, we have 


, , fll7oir/V||^||a: 

~ [ || 7 oir^^e'^'l^l 0 ||fe||e°'l^l otherwise, 

=: j{x,y). 


if X, P > 0 or X = P, 


(A39) 


b. Final steps 

Inserting the bounds in Eq. (A23) and Eq. (A39) into Eq. (A21) and letting c such that -^/HPj^X^ < c, we have 

AHt,L) 


< 


8\\Onh\\-Jo Jo 


dx Jdyj{x,x- y)(^- 


\Pl\\\Xl-i,l\ 


i+r1 




{cyf 


(2n)! /’ 


where 


f dx /'dy 7 (x,x - y)y" = IItoII^/^ [dx fdye'^J^ + ll^ll /dx [dye'^' 
Jo Jo Jo Jo Jo Jo 

l7oir/^+ 11^11 


i^-y) 


qCX _ I 


-y 


(A40) 


(A41) 


(n 0 l)(n 0 2) V. " " c' 

and we may take c' —> 0 if X, P > 0 or X = P. For L even (odd) we have [= P/2 (f^] = and [= Ll2 
= A^). Hence, for P = 1 the bound in the theorem follows. Finally, for P cx 1 the second term in Eq. (A40) vanishes 
and we have P = 0. 


Suppose 



(A 31 ) 


and let Q:( 0 ) = ao > 0 . Then 


a(x) 


Q:0 + 



— a(y) < ao + 

ay 


f dyf{y)s/a{y) =: /S{x). 

Jo 


(A 32 ) 


Lete > 0 and define 7 (ai) = l 3 (x) + e. Then 7 ( 21 ) > 7 ( 0 ) = ao + e > 0 , 
i.e., is differentiable and 

s/MJyjX = s/Axj = \/«o + e + f J‘y J V'yiy)^ (A 33 ) 

Jo di/ 

where 

= I j}y < \ f}yf(y)- (A34) 

As e > 0 was arbitrary, we hence have 

\/aix) < V^Q + 7 / dyf{y). (A 35 ) 

^ Jo 
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c. Correlation matrix 


The upper bound on the correlations a,b G {x,p}, in Eq. (A38) may be altered to allow for divergences at infinity by 

recalling that (see Eqs. (A25,A37)) 


lxx{y) lxp{y) \ _ p-crHuV^ (p-<rHBy\t 
7p.(2/) 7 pp( 2/) j ^ 


i.e. (see Eq. (A37)), 


Cxx (y) C^xp (y) ^ ^ ( ^xx (y) ^px {y)\ 
C-px (y) Cpp (y) / \ ^xpiy) ^pp(,y) J 


^xx{y) — ^xx {y)['yxx{0)cl.^{y) + 7xp(0)4p(t/)] + Ca;p{y)[ypx{0)cl.^iy) + 7pp(0)4p(y)], 
Ippiy) — ^px{,y')\')xx{^)Cpx{y) T7a;p(0)Cpp(?/)] + Cpp(y)[7p2;(0)Cp3.(y) + 7pp(0)Cpp(2/)]; 


(A42) 


(A43) 


and again using the bounds obtained in [16] which, however, increases the value of c in the bound. 


d. Multiple baths 


Eor some applications in quantum biology and condensed matter physics, one has a quantum system coupled to N baths 
which, using the Particle or Phonon mapping, can be written in the form 


N 


N 


^mul. 


(A44) 


where 


-| t 

A<"> = 5 E [ 


- (m)-r^(m) ^ (m) . ^(m) .-(n 

Xl/x) '+pI ’Plj'p) 


(A45) 


and G R, € [R- As in the rest of this work so far we assume that they couple only nearest 

^ij Ji^ ^iJ ^ 


neighbours, i.e. = 0 for \i — j\ > 1. We can truncate the N chains such that the to* chain contains modes; 


j(m) 


N 


N 


pmul. 


(A46) 


where 


im-l 


- 2 


E \Xm) , Am) p{m) Jm) 

[Xi ^i,j Xj -r Pi Pj 


ij=0 


E[A'{xu.,A'+t'r‘(PLj.,p\ 


(m) 

J 


(A47) 




where and Pp^ are principle submatrices of and P^™) corresponding to the non-truncated modes. These definitions 
are in analogy with those at the beginning of section A but generalised to the case of N non identical copies of the bath. 

Corollary 1 (Multiple chains) Lef pmui.^ pmui. above, such that UPl^AT^.,,, and 

max |jp(™)||} < c'j^. Then the error in truncating by is bounded by 


*] — tr[Oe 


Poe' 


N 

“*’*] I < Pi'n^,t,Lm) 


(A48) 


m—1 


A{L,t) ■= |tr[Oe 
where we have defined F > 0 as 

F^m,t,L) (A49) 

Cm. C,ry-, Cm. ^ \ J-J ~\~ L)\ \ Crwyi ' 


If X^'^\ p("*) > 0 or = p(™), we may take —>■ 0. If p('^') (x 1, we may replace ^7 ^‘^( 2 L +i)! • 


Jm) (m) 


7 o™^ = \Z)^ 44 )’ = tr[ 4 “^i 44 o], [ 7 p 4 ]i 7 = = tr[x\'^^p\"^^go], (A50) 

V 7ip 7pp / 
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collects the two-point m’^ bath correlations in the initial state of the whole system. Note that || < || and \ \P^'^ || < 

||PM||. 

As with theorem 2, one can allow for the two-point correlations collected in 7 g"*^ to diverge and still get a bound on F, see 
Section Ale. Often, one encounters systems interacting with multiple baths. We generalize to this setting in Section Aid. 
Proof. Starting from 


tr [Oi 


.-iff” 




*] -tr[d, 


.-iff? 




(A51) 


we add and subtract 

tr[de-‘*‘'-“*Poe*^‘™‘] (A52) 

A^ — 1 times where Htm. corresponds to but with some of the N baths truncated. Each time it is added and subtracted, 

different baths should be truncated. We then group the terms in pairs of 2 and redefine the system in each pair such that the 
system contains N — 1 baths (some truncated, some not). This step relies crucially on the fact that the system Hamiltonian Hg 
in not necessarily bounded. We then take the absolute value and apply the triangle inequality to the pairs followed by applying 
theorem 1 to each pair. q 

Thus the error introduced by truncating N chains, is bounded by the sum of the errors of truncating each chain individually. 
The explicit forms of the bound for the Particle and Phonon mapping can be found in section B. 


Appendix B: Derivation of the particle and phonon mapping chain trnneation bonnds 

From [9] we find Particle and Phonon mappings of Eq. (7) to be 

OO 

H =Hs + sj /3o(0)A5(6o(0) + (>q(0)) -f ^a„(0)6jj(0)&„(0) + sj /3„+i(0)(&Jj_|_]^(0)6„(0) + h.c.)^ (Bl) 

n—0 


and 


OO r1^ 1 

H =Hs + sj /3o(l)Asa;o(l) + ^ ^—^n(l) + + \/ /?ra+i(l)a;n(l)^n-i-i(l)) j (B2) 

n—0 

respectively. 6j,(0), (&„(0)), are creation (annihilation) operators. Define position and momentum operators for the particle 
mapping i„(0) := (6},(0) + 6„(0))/y2, p„(0) := i(&J^(0) - hn{Q))/V2 andi„(l) := '■= Pn^^)/s/^max 

for the phonon mapping. Write Eqs. (Bl), (B2) in terms of these new operators and compare these Eqs. with Eqs. (2) and (1). 
From here, together with the definition of the Jacobi matrices J{dX^) (see Eq. (162) in [9]), we find: 

For the particle mapping 


X = P = J{dX°), h= ^/2^o{0) As, dX°{x) = J{x)dx/7T. 


(B3) 


For the phonon mapping 


JjdXA 

^max 


P = 1 LCn 


h = 


//3o(l) 


UJr] 


As, 


dX^{x) = Jf\fx)dxl’K. 


From Eqs (15,156,160) in [9], 


m) 



dxJ{x)/Tt, 


r^max 

/3o(l) = / dxJis/x)/!:. 


(B4) 


(B5) 


Since the spectrum of a Jacobi matrix is equal to its minimally closed support interval [26], we have for the particle and phonon 
mappings: ||Ar|| = ||P|| = -^/H ArP|| = ccmax, and X > 0 iff oJmin > 0. For the Particle mapping we can use Eq. (A6) with 
c = P = cJmax to achieve 

A^(t,L) < + 7 +mI|A|I<) (B6) 
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where /tq given by Eq. (11). For the massive Phonon mapping, we replace 8 with 4, po with /ri, and /(^ + 1)! by 

+ 1)! in Eq. (B 6 ). For the massless Phonon chain mapping, we use Eq. (A 6 ) again, to achieve 


^max y.tiL/-\- ij! 


1) 


1/2 


- 1 


(B7) 


We can write the 70 matrix for the Phonon mapping in terms of the original Xn and pn coordinates of Eq. 10, to find 


7o 


^max jxx 


"ixp 


Ipx 


:7pp 


[yab]n,i = tr[a„6iPo]- 


(B 8 ) 


Appendix C: Fock space truncation 

In this section we derive bounds on the error introduced by truncating the local Hilbert spaces of the bath. To this end, we 
define the projector 


1™. — Imo ® Imi ® ■ ■ ■ I (Cl) 

where acts on the z’th site of the bath and truncates the local Hilbert space according to 

771 

= (C 2 ) 

n =0 


For bounded observables acting on the system O, ||0|| < 00 , we consider 

A^(f) = |tr[de-'‘*poe'‘*] -tr[6e-*‘^-poe‘‘^"‘] 
i.e., the error introduced by evolving the system according to 

Hm = tmiltm = Hs+HB+Vm 


(C3) 


(C4) 


instead of H = Hs + Hb + V. Here, with the notation = trnXitm and p^ = ImPitm, the individual terms read 
Vrn = ^ < 8 > ^ 0 ™ and 


= 1. 






(C5) 


where we note that and ^ while for i ^ j we do have lm.XiXjlm. = xY^xJ^ and 

ImPiPjffm = pTpT- denote 


U{t) = 


(C 6 ) 


Proceeding as in Eqs. (A9-A12), we find 

^m{t) = |tr[de"'*^poe**^] - tr[6e"**^"‘poe“^”*]| < 2||d|l-\/tr[[[7t(f) - Uln{t)][U{t) - Um{t)]Qo\, (C7) 


where now, as U^{t)Um{t) = 15 ® 1 ^, 




tr[(l - lrn) 0 o] - 2 3? J da; — tr[[/'l'(x)(7r„(x)po] 


(C 8 ) 


and 


-i^U\x)Urr.ix) = 
dx 


-ixHn _ 




— ixHrt 


(C9) 
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FIG. 3. Fock space truncation error (Eq. (17)) for the model in Eq. (19) using the particle mapping with A/cuc = 1, a = 0.8, s = 1 for initial 
state ^0 = ® = I t)(t I ™d the vacuum. We tmncate each local Flilbert space at the same value nii = m and L has the values 

3 to 6, but are indistinguishable (e.g. the difference between the L = 6 and L = 3 curve at the point denoted by a square is 4.86 x 10“®). 
Lines are guides to the eye. 


The Cauchy-Schwarz inequality yields 


^ — tT[U\x)Umix)go]f<tT[h'^e =: e^(x) 


such that 


A^(f) <4||df ^tr[(l - l^)eo] +‘2. j dx^/e^x)^ . 

The error em(x) may be obtained numerically: We have 

= 6 ““^™ W{x)tmW{x)e^^^^ + w(x){1 - 


(CIO) 


(Cll) 


(C12) 


such that, recalling Eq. (A25), i.e., that xo{t) = e'^^’^xoe = X]/c ^o\{^)Pk, the computation of em[x) is 


— IxHrn A r-X^Hn 


reduced to obtaining the coefficients CQ^(f) and CQ^^(t) and expectations in e observables of the form 


and 


for r, s G {x,p}. 


- xJP) 




(C13) 


(C14) 


Appendix D: Further numerical examples of Fock space truncation 

In this section, we give the chain coefficients used in the numerical simulations for the Fock space truncation and analyse 
further the numerical results. We start with the particle mapping of the spin-boson model: From Eq. (B3), we have the relation 
between the X and P matrices and the Jacobi matrix. The coefficients of the Jacobi matrix for the spin-boson spectral density 
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of Eq. (19) can be found in [7] or [9]. From Eqs. (248), (249) in [9], we have for n S N°, 

X f, I ^ 

2 [ (s + 2n)(2 + s + 2n)J ’ 

2fn,n+l — ^n+l,n 

a;c(l + n)(l + s + n) /S + s + 2 n 
(5 4-2-1- 27r)(3 4- 5 4- 27r) V 1 4- 5 4- 27z 

and all other matrix elements zero. /3o(0) can be found in Eq. (250) in [9], thus from Eq. (B3), we have 


(Dl) 

(D2) 

(D3) 


h = 


UJc 



(D4) 


We can now do the same for the phonon mapping written in terms of and p . Using Eq. (B4), we have the relation between 
the X and P matrices and the Jacobi matrix. From [9], we obtain the coefficients of the Jacobi matrix. Thus, using Eqs. (255), 
(256) in [9], we have for n G N°, 


^n,n+l ^n+l,n 


Wc2(l- 

JTTT- 


-n)(Z 


4n)(6 4- s 4- 4n) \ 2 + s 


(D5) 

(D6) 


and all other matrix elements zero. /3o(l) can be found in Eq. (257) in [9], thus from Eq. (B4), we have 


h 


= UJc 



(D7) 


The results for the particle mapping with Ohmic spectral density are plotted in Fig. 3 and for the phonon mapping in Fig. 2. In 
both cases, the plots suggest that the super ohmic spectral densities have smaller truncation error. For the particle mappings, 
in the plots we have probed the zero Kelvin thermal state (which corresponds to the chain vacuum state), where as for the 
phonon mappings we have probed a squeezed vacuum state, which is highly populated. We see that the error has slightly worse 
decay with increasing m than in the particle mappings cases. This is intuitively what one would expect, since more of the bath 
population is being truncated. 


We probed the vacuum state of the chain, which corresponds to a squeezed 
vacuum state of the continuous bath of hai'monic oscillators. This will be 


shown in an upcoming article [27]. 















